{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.7 Nonlinear problems\n",
    "We want to solve a nonlinear PDE. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A simple scalar PDE\n",
    "We consider the simple PDE \n",
    "\n",
    "$$\n",
    "- \\Delta u + 3 u^3 = 1 \\text{ in } \\Omega\n",
    "$$\n",
    "\n",
    "on the unit square $\\Omega = (0,1)^2$. \n",
    "\n",
    "We note that this PDE can also be formulated as a nonlinear minimization problem (cf. [3.8](../unit-3.8-nonlmin/nonlmin.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from netgen import gui\n",
    "from netgen.geom2d import unit_square\n",
    "from ngsolve import *\n",
    "\n",
    "mesh = Mesh (unit_square.GenerateMesh(maxh=0.3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In NGSolve we can solve the PDE conveniently using the *linearization* feature of `SymbolicBFI`.\n",
    "\n",
    "The `BilinearForm` (which is not bilinear!) needed in the weak formulation is\n",
    "$$\n",
    "  A(u,v) = \\int_{\\Omega} \\nabla u \\nabla v + 3 u^3 v - 1 v ~ dx \\quad ( = 0 ~ \\forall~v \\in H^1_0)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "V = H1(mesh, order=3, dirichlet=[1,2,3,4])\n",
    "u,v = V.TnT()\n",
    "a = BilinearForm(V)\n",
    "a += (grad(u) * grad(v) + 3*u**3*v- 1 * v)*dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Newton's method\n",
    "\n",
    "We use Newton's method and make the loop:\n",
    "\n",
    "* Given an initial guess $u^0$\n",
    "* loop over $i=0,..$ until convergence:\n",
    "  * Compute linearization: $A u^i + \\delta A(u^i) \\Delta u^{i} = 0$:\n",
    "    * $f^i = A u^i$ \n",
    "    * $B^i = \\delta A(u^i)$ \n",
    "    * Solve $B^i \\Delta u^i = -f^i$\n",
    "  * Update $u^{i+1} = u^i + \\Delta u^{i}$\n",
    "  * Evaluate stopping criteria\n",
    "\n",
    "As a stopping criteria we take $\\langle A u^i,\\Delta u^i \\rangle = \\langle A u^i, A u^i \\rangle_{(B^i)^{-1}}< \\varepsilon$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def SimpleNewtonSolve(gfu,a,tol=1e-13,maxits=25):\n",
    "    res = gfu.vec.CreateVector()\n",
    "    du = gfu.vec.CreateVector()\n",
    "    fes = gfu.space\n",
    "    for it in range(maxits):\n",
    "        print (\"Iteration {:3}  \".format(it),end=\"\")\n",
    "        a.Apply(gfu.vec, res)\n",
    "        a.AssembleLinearization(gfu.vec)\n",
    "        du.data = a.mat.Inverse(fes.FreeDofs()) * res\n",
    "        gfu.vec.data -= du\n",
    "\n",
    "        #stopping criteria\n",
    "        stopcritval = sqrt(abs(InnerProduct(du,res)))\n",
    "        print (\"<A u\",it,\", A u\",it,\">_{-1}^0.5 = \", stopcritval)\n",
    "        if stopcritval < tol:\n",
    "            break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(V)\n",
    "Draw(gfu,mesh,\"u\")\n",
    "SimpleNewtonSolve(gfu,a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "There are also some solvers shipped with NGSolve now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve.solvers import *\n",
    "help(Newton)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "gfu.vec[:]=0\n",
    "Newton(a,gfu,freedofs=gfu.space.FreeDofs(),maxit=100,maxerr=1e-11,inverse=\"umfpack\",dampfactor=1,printing=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A trivial problem:\n",
    "\n",
    "$$\n",
    "  5 u^2 = 1, \\qquad u \\in \\mathbb{R}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "V = NumberSpace(mesh)\n",
    "u,v = V.TnT()\n",
    "a = BilinearForm(V)\n",
    "a += ( 5*u*u*v - 1 * v)*dx\n",
    "gfu = GridFunction(V)\n",
    "gfu.vec[:] = 1\n",
    "SimpleNewtonSolve(gfu,a)\n",
    "\n",
    "print(\"\\nscalar solution\", gfu.vec[0], \"(exact: \", sqrt(0.2), \")\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Another example: Stationary Navier-Stokes:\n",
    "Find $\\mathbf{u} \\in \\mathbf{V}$, $p \\in Q$, $\\lambda \\in \\mathbb{R}$ so that\n",
    "\\begin{align}\n",
    "\\int_{\\Omega} \\nu \\nabla \\mathbf{u} : \\nabla \\mathbf{v} + (\\mathbf{u} \\cdot \\nabla) \\mathbf{u} \\cdot \\mathbf{v}& - \\int_{\\Omega} \\operatorname{div}(\\mathbf{v}) p & &= \\int \\mathbf{f}  \\cdot \\mathbf{v}  && \\forall \\mathbf{v} \\in \\mathbf{V}, \\\\ \n",
    "- \\int_{\\Omega} \\operatorname{div}(\\mathbf{u}) q & & \n",
    "+ \\int_{\\Omega} \\lambda q\n",
    "&= 0 && \\forall q \\in Q, \\\\\n",
    "& \\int_{\\Omega} \\mu p & &= 0 && \\forall \\mu \\in \\mathbb{R}.\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mesh = Mesh (unit_square.GenerateMesh(maxh=0.05)); nu = Parameter(1)\n",
    "V = VectorH1(mesh,order=3,dirichlet=\"bottom|right|top|left\")\n",
    "Q = H1(mesh,order=2); \n",
    "N = NumberSpace(mesh); \n",
    "X = FESpace([V,Q,N])\n",
    "(u,p,lam), (v,q,mu) = X.TnT()\n",
    "a = BilinearForm(X)\n",
    "a += (nu*InnerProduct(grad(u),grad(v))+InnerProduct(grad(u)*u,v)\n",
    "      -div(u)*q-div(v)*p-lam*q-mu*p)*dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(X)\n",
    "gfu.components[0].Set(CoefficientFunction((4*x*(1-x),0)),\n",
    "                      definedon=mesh.Boundaries(\"top\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "SimpleNewtonSolve(gfu,a)\n",
    "Draw(gfu.components[1],mesh,\"p\")\n",
    "Draw(gfu.components[0],mesh,\"u\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "nu.Set(0.01)\n",
    "SimpleNewtonSolve(gfu,a)\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "nu.Set(0.001)\n",
    "SimpleNewtonSolve(gfu,a)\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "nu.Set(0.001)\n",
    "gfu.components[0].Set(CoefficientFunction((4*x*(1-x),0)),definedon=mesh.Boundaries(\"top\"))\n",
    "Newton(a,gfu,maxit=20,dampfactor=0.1)\n",
    "Redraw()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
